!--------------------------------
!--------------------------------
!--------------------------------
subroutine grid_uniform
  use rbf_type
  
  use variable, only:xmax,ymax,zmax,ddx,ddy,ddz,ddx2,ddy2,ddz2,ry,rz,beta,nip,iddx2,iddy2,iddz2
  use variable, only:Nx,Ny,Nz,dddx,dddy,dddz,nn,wall,plan
  use variable, only: confine,surfactant
  use variable, only:inside,nee_cube,xm,nematic_l
  use variable, only:surface,turns,width,Lw,Gw,Mw
  use variable, only:idxdy,idxdz,idydz
  implicit none
 
  real(dp)      :: x,y,z,r,checkx,checky,checkz,check1,check2,dplan,xmm,radiusp,ymm
  real(dp)      :: lstri,rstri,lgap,rgap
  real(dp), dimension(3) :: xs
  integer (i4b)  :: i,j,k,cont,cont_nn,cont3,cont_total,pl
  real(dp), dimension(:),allocatable :: xm_grid,xm_nn
  logical , dimension(:),allocatable :: plan_tmp,wall_tmp


  write(*,*)'------------------------------------------------------------'
  xmax =dble(beta*2.d0)         ! define the lenght of the simulation box
  ymax =dble(beta*2.d0)*ry 
  zmax =dble(beta*2.d0)*rz 

  write(*,160)"   xmax   =",xmax , '   Lenght x (nm)  =',xmax*nematic_l
  write(*,160)"   ymax   =",ymax , '   Lenght y (nm)  =',ymax*nematic_l
  write(*,160)"   zmax   =",zmax , '   Lenght z (nm)  =',zmax*nematic_l
  write(*,*)'------------------------------------------------------------'

  ddx = xmax/dble(Nx-1)         ! Distance between points
  ddy = ymax/dble(Ny-1)
  ddz = zmax/dble(Nz-1)

  write(*,150)"   Mesh   ddx  =",ddx
  write(*,150)"   Mesh   ddy  =",ddy
  write(*,150)"   Mesh   ddz  =",ddz
  write(*,*)'------------------------------------------------------------'

150  format (a,f10.4)
160  format (a,f12.4,a,f12.4)

  iddx2 =  1.d0/ddx**2
  iddy2 =  1.d0/ddy**2
  iddz2 =  1.d0/ddz**2

  ddx2 =  ddx**2
  ddy2 =  ddy**2
  ddz2 =  ddz**2

  idxdy = 1.d0/(4.d0*ddx*ddy)
  idxdz = 1.d0/(4.d0*ddx*ddz)
  idydz = 1.d0/(4.d0*ddy*ddz)

  dddx = 0.5d0/ddx 
  dddy = 0.5d0/ddy
  dddz = 0.5d0/ddz

  xs(1) = 0.5d0*xmax
  xs(2) = 0.5d0*ymax
  xs(3) = 0.5d0*zmax

  R = beta*0.5d0

!******************************************************
!loop for surface nodes before internal nodes loop
  allocate(inside   (Nx,Ny,Nz))
  allocate(nee_cube (Nx,Ny,Nz))
  allocate(xm_grid  (nx*ny*nz*3))
  allocate(surface  (Nx,Ny,Nz))
  allocate(xm_nn    (nx*ny*nz*3))
  allocate(wall_tmp(Nx*Ny*Nz))
  allocate(plan_tmp(Nx*ny*nz))

  nee_cube = 0
  surface=.false.
  inside=.false.
  cont       = 0
  cont_nn    = 0
  cont_total = 0  
  cont3      = 0
  width      = Mw/nematic_l
  dplan = (xmax - dble(turns)*width)/dble(turns)
     
  
  write(*,*)'dplan  =',dplan
  write(*,*)'width  =',width
  write(*,*)
  write(*,*)'  Middle line (in nm) = ', Mw 
  write(*,*)'  Gap  Space  (in nm) = ', Gw
  write(*,*)'  Line  width (in nm) = ', Lw

  Mw = Mw /nematic_l
  Gw = Gw /nematic_l
  Lw = Lw /nematic_l
  write(*,*)
  write(*,*)'  Middle line  = ', Mw 
  write(*,*)'  Gap  Space   = ', Gw
  write(*,*)'  Line  width  = ', Lw

!---------------------------------------------------------------
     do i=1,Nx
        do j=1,Ny
           do k=1,Nz

              x = ddx*(i-1)
              y = ddy*(j-1)
              z = ddz*(k-1)
             
             
              if  (k.eq.1.or.k.eq.NZ)then
                 cont_nn = cont_nn +1 
              
                 cont3 = cont3 + 1  
                 surface(i,j,k)=.true.
                 inside(i,j,k)=.false.
             
                 xm_nn(3*cont_nn - 2) =  x - xs(1)
                 xm_nn(3*cont_nn - 1) =  y - xs(2)
                 xm_nn(3*cont_nn    ) =  z - xs(3)
  
                 nee_cube(i,j,k)=cont3
              
                 if(surfactant.eq.2)then  ! Stripes
                    if (k.eq.1)then
                       xmm = x 
                       do pl=0,turns-1
                          if (xmm.gt. (pl*(width+dplan)+dplan*0.5) &
                         .and.xmm.lt.((pl*(width+dplan)+dplan*0.5) + width))then
                             plan_tmp(cont_nn)=.true.
                          endif
                       enddo
                    endif
                 else if(surfactant.eq.5)then   ! Chevron

                 else if(surfactant.eq.6) then  ! T-shape

                 else if(surfactant.eq.7)then   ! Rectangles
                    
                 else if(surfactant.eq.8)then   ! Crossroads
                    if (k.eq.1)then
                       xmm = x 
                       do pl=0,turns-1
                          if (xmm.gt. (pl*(width+dplan)+dplan*0.5).and.xmm.lt.((pl*(width+dplan)+dplan*0.5) + width))then
                             plan_tmp(cont_nn)=.true.
                          endif
                       enddo
                    endif

                    if (k.eq.1)then
                       xmm = x
                       ymm = y
                       do pl=0,turns-1
                          lgap = ((pl*(width+dplan)+dplan*0.5) - Gw)
                          rgap =((pl* (width+dplan)+dplan*0.5) + width+  Gw)
                          lstri =(  pl*(width+dplan)+dplan*0.5)
                          rstri = ((pl*(width+dplan)+dplan*0.5) + width)
                          
                          if ((-ymm+xmm).le.(Lw*0.5) .and.(-ymm+xmm).ge.(-Lw*0.5)  )then
                             plan_tmp(cont_nn)=.true.      
                          endif
                          if ((xmm.gt.lgap.and.xmm.le.lstri).or.(xmm.lt.rgap.and.xmm.ge.rstri))then  
                             plan_tmp(cont_nn)=.false.      
                          endif
                          
                       enddo
                    endif
                    

                 else if(surfactant.eq.4)then   ! Circles
                    if (k.eq.1)then
                       radiusp = xm_nn(3*cont_nn - 2)**2 + xm_nn(3*cont_nn - 1)**2
                       if (radiusp.gt.width**2.and.radiusp.lt.(width+10.d0)**2)then
                          plan_tmp(cont_nn)=.true.
                       endif
                    endif
                 
                 endif

           endif
        enddo
     enddo
  enddo

  !********************************************

  do i=1,Nx
     do j=1,Ny
        do k=1,Nz
        
           x = ddx*(i-1)
           y = ddy*(j-1)
           z = ddz*(k-1)
           
           if (.not.surface(i,j,k))then
              
              cont3 = cont3 + 1
              cont = cont + 1

              xm_grid(3*cont - 2) =  x - xs(1)
              xm_grid(3*cont - 1) =  y - xs(2)
              xm_grid(3*cont    ) =  z - xs(3)
              inside(i,j,k)=.true.
              nee_cube(i,j,k)=cont3
              if (k.eq.2) then
                 wall_tmp (cont) = .true.
              endif

              if (k.eq.(Nz-1)) then
                 wall_tmp(cont) =.true.
              endif
           endif

        enddo
     enddo
  enddo
!********************************************


147 format (3I5,3f12.6)
  
 
 nip = cont
 nn = cont_nn

!print *, 'nn+nip=', nn+nip, 'cont3=', cont3


 allocate (xm(nn*3 + nip*3))
 allocate (plan(nn))
 allocate (wall(nip+nn))
 xm(1 : 3*nn )                  = xm_nn (1:3*cont_nn)
 xm((3*nn)+1 : (3*nn)+(3*nip) ) = xm_grid(1:3*cont)
 plan(1:nn) = plan_tmp(1:nn)
 wall(nn+1:nn+nip)  = wall_tmp (1:cont)
 deallocate (xm_nn,xm_grid,plan_tmp,wall_tmp)


end subroutine grid_uniform
!--------------------------------------------------------------------
!--------------------------------------------------------------------
!--------------------------------------------------------------------
subroutine cube
  use rbf_type
 use variable, only:inside,nee_cube,nx,ny,nz,surface,nodes
  use variable, only:ix1_surf,ix2_surf
 use variable, only: ix1,ix2,ix3,ix4,ix5

  implicit none

!-----------------------------------------------------------------
  integer(i4b)  :: i, j, k,i0
  integer (i4b) :: ii,jj,kk,NN
  
 allocate(ix1      (6,nodes))
 allocate(ix2      (6,nodes))
 allocate(ix3      (6,nodes))
 allocate(ix4      (6,nodes))
 allocate(ix5      (6,nodes))

 allocate(ix1_surf (6,nodes))
 allocate(ix2_surf (6,nodes))

 ix1=0
 ix1_surf=0
 ix2_surf=0

  do i=1, nx
     do j=1, ny
        do k=1, nz

              i0=nee_cube(i ,j ,k )
              
              if(i.ge.2)then
                 
                 ii = i-1
                 NN = 1
                               ix1  (NN,i0)=nee_cube(ii,j  ,k  )

                 if(j.ge.2   ) ix2  (NN,i0)=nee_cube(ii,j-1,k  )
                 if(j.eq.1   ) ix2  (NN,i0)=nee_cube(ii,Ny ,k  )

                 if(j.le.Ny-1) ix3  (NN,i0)=nee_cube(ii,j+1,k  )
                 if(j.eq.Ny  ) ix3  (NN,i0)=nee_cube(ii,1  ,k  )

                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(ii,j  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(ii,j  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(ii,j  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(ii,j  ,1  )
    
              else if(i.eq.1)then

                 ii = Nx
                 NN = 1
                               ix1  (NN,i0)=nee_cube(ii ,j  ,k  )

                 if(j.ge.2   ) ix2  (NN,i0)=nee_cube(ii ,j-1,k  )
                 if(j.eq.1   ) ix2  (NN,i0)=nee_cube(ii ,Ny ,k  )

                 if(j.le.Ny-1) ix3  (NN,i0)=nee_cube(ii ,j+1,k  )
                 if(j.eq.Ny  ) ix3  (NN,i0)=nee_cube(ii ,1  ,k  )
                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(ii ,j  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(ii ,j  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(ii ,j  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(ii ,j  ,1  )

              endif

!--------------------------------------------------------------------------
!--------------------------------------------------------------------------

              if(i.le.Nx-1)then
                 
                 ii = i+1
                 NN = 2
                               ix1  (NN,i0)=nee_cube(ii,j  ,k  )

                 if(j.ge.2   ) ix2  (NN,i0)=nee_cube(ii,j-1,k  )
                 if(j.eq.1   ) ix2  (NN,i0)=nee_cube(ii,Ny ,k  )

                 if(j.le.Ny-1) ix3  (NN,i0)=nee_cube(ii,j+1,k  )
                 if(j.eq.Ny  ) ix3  (NN,i0)=nee_cube(ii,1  ,k  )

                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(ii,j  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(ii,j  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(ii,j  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(ii,j  ,1  )
    
              else if(i.eq.Nx)then

                 ii = 1
                 NN = 2
                               ix1  (NN,i0)=nee_cube(ii ,j  ,k  )

                 if(j.ge.2   ) ix2  (NN,i0)=nee_cube(ii ,j-1,k  )
                 if(j.eq.1   ) ix2  (NN,i0)=nee_cube(ii ,Ny ,k  )

                 if(j.le.Ny-1) ix3  (NN,i0)=nee_cube(ii ,j+1,k  )
                 if(j.eq.Ny  ) ix3  (NN,i0)=nee_cube(ii ,1  ,k  )
                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(ii ,j  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(ii ,j  ,Nz )
                                  
                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(ii ,j  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(ii ,j  ,1  )

              endif


!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!---------------------------------------------------------------------------

              if(j.ge.2)then
                 
                 jj = j-1
                 NN = 3
                               ix1  (NN,i0)=nee_cube(i  ,jj  ,k  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1,jj  ,k  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx ,jj  ,k  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1,jj  ,k  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1  ,jj  ,k  )

                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(i  ,jj  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(i  ,jj  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(i  ,jj  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(i  ,jj  ,1  )
    
              else if(j.eq.1)then

                 jj = Ny
                 NN = 3
                               ix1  (NN,i0)=nee_cube(i   ,jj  ,k  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1 ,jj  ,k  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx  ,jj  ,k  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1 ,jj  ,k  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1   ,jj  ,k  )
                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(i   ,jj  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(i   ,jj  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(i   ,jj  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(i   ,jj  ,1  )

              endif

!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!---------------------------------------------------------------------------

              if(j.le.Ny-1)then
                 
                 jj = j+1
                 NN = 4
                               ix1  (NN,i0)=nee_cube(i  ,jj  ,k  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1,jj  ,k  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx ,jj  ,k  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1,jj  ,k  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1  ,jj  ,k  )

                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(i  ,jj  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(i  ,jj  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(i  ,jj  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(i  ,jj  ,1  )
    
              else if(j.eq.Ny)then

                 jj = 1
                 NN = 4
                               ix1  (NN,i0)=nee_cube(i   ,jj  ,k  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1 ,jj  ,k  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx  ,jj  ,k  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1 ,jj  ,k  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1   ,jj  ,k  )
                 
                 if(k.ge.2   ) ix4  (NN,i0)=nee_cube(i   ,jj  ,k-1)
                 if(k.eq.1   ) ix4  (NN,i0)=nee_cube(i   ,jj  ,Nz )
                                  

                 if(k.le.Nz-1) ix5  (NN,i0)=nee_cube(i   ,jj  ,k+1)
                 if(k.eq.Nz  ) ix5  (NN,i0)=nee_cube(i   ,jj  ,1  )

              endif


!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------

              if(k.ge.2)then
                 
                 kk = k-1
                 NN = 5
                               ix1  (NN,i0)=nee_cube(i   ,j    ,kk  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1 ,j    ,kk  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx  ,j    ,kk  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1 ,j    ,kk  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1   ,j    ,kk  )

                 if(j.ge.2   ) ix4  (NN,i0)=nee_cube(i   ,j-1  ,kk  )
                 if(j.eq.1   ) ix4  (NN,i0)=nee_cube(i   ,Ny   ,kk  )

                 if(j.le.Ny-1) ix5  (NN,i0)=nee_cube(i   ,j+1  ,kk  )
                 if(j.eq.Ny  ) ix5  (NN,i0)=nee_cube(i   ,1    ,kk  )

    
              else if(k.eq.1)then

                 kk = Nz
                 NN = 5

         !                      ix1  (NN,i0)=nee_cube(i  ,j   ,kk  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1,j   ,kk  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx ,j   ,kk  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1,j   ,kk  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1  ,j   ,kk  )

                 if(j.ge.2   ) ix4  (NN,i0)=nee_cube(i  ,j-1 ,kk  )
                 if(j.eq.1   ) ix4  (NN,i0)=nee_cube(i  ,Ny  ,kk  )

                 if(j.le.Ny-1) ix5  (NN,i0)=nee_cube(i  ,j+1 ,kk  )
                 if(j.eq.Ny  ) ix5  (NN,i0)=nee_cube(i  ,1   ,kk  )

              endif

!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------
!-----------------------------------------------------------------------------------

              if(k.le.Nz-1)then
                 
                 kk = k+1
                 NN = 6
                               ix1  (NN,i0)=nee_cube(i  ,j   ,kk  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1,j   ,kk  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx ,j   ,kk  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1,j   ,kk  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1  ,j   ,kk  )

                 if(j.ge.2   ) ix4  (NN,i0)=nee_cube(i  ,j-1  ,kk  )
                 if(j.eq.1   ) ix4  (NN,i0)=nee_cube(i  ,Ny   ,kk  )

                 if(j.le.Ny-1) ix5  (NN,i0)=nee_cube(i  ,j+1  ,kk  )
                 if(j.eq.Ny  ) ix5  (NN,i0)=nee_cube(i  ,1    ,kk  )

    
              else if(k.eq.Nz)then

                 kk = Nz
                 NN = 6

          !                     ix1  (NN,i0)=nee_cube(i  ,j   ,kk  )

                 if(i.ge.2   ) ix2  (NN,i0)=nee_cube(i-1,j   ,kk  )
                 if(i.eq.1   ) ix2  (NN,i0)=nee_cube(Nx ,j   ,kk  )

                 if(i.le.Nx-1) ix3  (NN,i0)=nee_cube(i+1,j   ,kk  )
                 if(i.eq.Nx  ) ix3  (NN,i0)=nee_cube(1  ,j   ,kk  )

                 if(j.ge.2   ) ix4  (NN,i0)=nee_cube(i  ,j-1 ,kk  )
                 if(j.eq.1   ) ix4  (NN,i0)=nee_cube(i  ,Ny  ,kk  )

                 if(j.le.Ny-1) ix5  (NN,i0)=nee_cube(i  ,j+1 ,kk  )
                 if(j.eq.Ny  ) ix5  (NN,i0)=nee_cube(i  ,1   ,kk  )

              endif




!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------


             if (surface(i,j,k))then

              i0=nee_cube(i ,j ,k )
              if(k.eq.1.or.k.eq.Nz)then
                 
                 if(i.gt.2)then
                    if( nee_cube(i-1,j,k).ne.0)     ix1_surf(1,i0)=nee_cube(i-1,j,k)
                    if( nee_cube(i-2,j,k).ne.0)     ix2_surf(1,i0)=nee_cube(i-2,j,k)
                 else if(i.eq.2)then
                    if( nee_cube(1  ,j,k).ne.0)     ix1_surf(1,i0)=nee_cube(1  ,j,k) 
                    if( nee_cube(Nx ,j,k).ne.0)     ix2_surf(1,i0)=nee_cube(Nx ,j,k)
                 else if(i.eq.1)then
                    if( nee_cube(Nx,j,k).ne.0)      ix1_surf(1,i0)=nee_cube(Nx ,j,k) 
                    if( nee_cube(Nx-1,j,k).ne.0)    ix2_surf(1,i0)=nee_cube(Nx-1,j,k)
                 endif

                 if(i.lt.(Nx-1))then
                    if( nee_cube(i+1,j,k).ne.0)     ix1_surf(2,i0)=nee_cube(i+1,j,k)
                    if( nee_cube(i+2,j,k).ne.0)     ix2_surf(2,i0)=nee_cube(i+2,j,k) 
                 else if(i.eq.Nx-1)then
                    if( nee_cube(Nx ,j,k).ne.0)     ix1_surf(2,i0)=nee_cube(Nx ,j,k)
                    if( nee_cube(1  ,j,k).ne.0)     ix2_surf(2,i0)=nee_cube(1  ,j,k)
                 else if(i.eq.Nx)then
                    if( nee_cube(1  ,j,k).ne.0)     ix1_surf(2,i0)=nee_cube(1,j,k)
                    if( nee_cube(2  ,j,k).ne.0)     ix2_surf(2,i0)=nee_cube(2,j,k)
                 endif


                 if(j.gt.2)then
                    if( nee_cube(i,j-1,k).ne.0)     ix1_surf(3,i0)=nee_cube(i,j-1,k)
                    if( nee_cube(i,j-2,k).ne.0)     ix2_surf(3,i0)=nee_cube(i,j-2,k)
                 else if(j.eq.2)then
                    if( nee_cube(i,1  ,k).ne.0)     ix1_surf(3,i0)=nee_cube(i,1  ,k)
                    if( nee_cube(i,Ny ,k).ne.0)     ix2_surf(3,i0)=nee_cube(i,Ny ,k)
                 else if(j.eq.1)then
                    if( nee_cube(i,Ny ,k).ne.0)     ix1_surf(3,i0)=nee_cube(i,Ny ,k)
                    if( nee_cube(i,Ny-1,k).ne.0)    ix2_surf(3,i0)=nee_cube(i,Ny-1,k)
                 endif


                 if(j.lt.(Ny-1))then
                    if( nee_cube(i,j+1,k).ne.0)     ix1_surf(4,i0)=nee_cube(i,j+1,k)
                    if( nee_cube(i,j+2,k).ne.0)     ix2_surf(4,i0)=nee_cube(i,j+2,k)
                 else if(j.eq.Ny-1)then
                    if( nee_cube(i,Ny ,k).ne.0)     ix1_surf(4,i0)=nee_cube(i,Ny ,k)
                    if( nee_cube(i,1  ,k).ne.0)     ix2_surf(4,i0)=nee_cube(i,1  ,k)
                 else if(j.eq.Ny)then
                    if( nee_cube(i,1  ,k).ne.0)     ix1_surf(4,i0)=nee_cube(i,1  ,k)
                    if( nee_cube(i,2  ,k).ne.0)     ix2_surf(4,i0)=nee_cube(i,2  ,k)
                 endif

                 if(k.eq.Nz)then
                   
                    if( nee_cube(i,j,k-1).ne.0)     ix1_surf(5,i0)=nee_cube(i,j,k-1) 
                    if( nee_cube(i,j,k-2).ne.0)     ix2_surf(5,i0)=nee_cube(i,j,k-2)                    
                 endif

                 if(k.eq.1)then

                    if( nee_cube(i,j,k+1).ne.0)     ix1_surf(6,i0)=nee_cube(i,j,k+1) 
                    if( nee_cube(i,j,k+2).ne.0)     ix2_surf(6,i0)=nee_cube(i,j,k+2) 

                 endif
                 
              endif
           endif


        end do
     end do
  end do


          
end subroutine cube
!-----------------------------------------------------
!-----------------------------------------------------
!-----------------------------------------------------
subroutine drop_normal
  use rbf_type
  use variable, only: nn, xm, normal,planar,beta,pi,tangential,tangential_2,ix1_surf
  implicit none
  integer(i4b) :: i, j, ii, jj
  real(dp)     :: r,zz,rr,phi,theta
  real(dp), dimension(3):: dx,dx_norm

  normal = 0.d0
 
  do i = 1, nn
     ii = 3*i
     dx = xm(ii - 2:ii)
     r  = sqrt( dot_product(dx,dx) )
     normal(ii - 2:ii) = dx/r

     
  end do

  do i = 1, nn
     ii = 3*i

    if (ix1_surf(6,i).ne.0) then         ! Z  + 1
     normal (ii-2) = 0.d0
     normal (ii-1) = 0.d0
     normal (ii  ) = 1.d0
  !     coef =  1.d0
  !     top =  .false.
    else if (ix1_surf(5,i).ne.0) then    ! Z - 1


!       coef = -1.d0
!       top = .true.
     normal (ii-2) = 0.d0
     normal (ii-1) = 0.d0
     normal (ii  ) = - 1.d0
    endif



 !    normal (ii-2) = 0.d0
 !    normal (ii-1) = 0.d0
 !    normal (ii  ) = 1.d0

     tangential (ii-2 )  = 1.0d0
     tangential (ii-1 )  = 0.0d0
     tangential (ii   )  = 0.0d0

     tangential_2 (ii-2 )=0.d0
     tangential_2 (ii-1 )=1.d0
     tangential_2 (ii   )=0.d0
  enddo
 

end subroutine drop_normal
